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Abstract 

Throughout  the  Air  Force  and  DoD,  chlorinated  aliphatic  hydrocarbons  (CAHs), 
which  are  a  type  of  Dense  Non- Aqueous  Phase  Liquid  or  “DNAPL,”  have  been 
frequently  used.  Their  disposal  has  led  to  subsurface  soil  and  groundwater 
contamination.  DNAPLs  move  through  soils  and  groundwater,  leaving  behind  residual 
separate  phase  contamination  as  well  as  pools  atop  low  penneability  layers.  DNAPL 
contamination  is  difficult  to  clean  up  with  existing  environmental  remediation 
technologies.  In  this  study  a  numerical  model  constructed  using  the  DoD  Groundwater 
Modeling  System  (GMS)  is  used  as  a  tool  to  assess  how  cracks  in  a  low  permeability 
layer,  either  pre-existing  or  resulting  from  the  interaction  between  a  DNAPL  pool  and  the 
layer,  might  ultimately  impact  the  longevity  of  a  dissolved  phase  CAH  plume  that  will  be 
generated  by  the  DNAPL  source.  The  conceptual  model  posits  a  DNAPL  source  in  a 
high  permeability  sand  layer  is  sitting  atop  a  low  permeability  clay  layer.  The  model 
simulates  DNAPL  dissolving  into  groundwater  flowing  through  the  sand,  and  assumes 
the  dissolved  CAH  is  transported  by  the  processes  of  advection,  dispersion,  and  sorption. 
These  transport  processes  in  the  sand  are  coupled  to  processes  occurring  in  the  low 
penneability  clay.  Dissolved  phase  CAH  is  assumed  to  be  transported  by  diffusion  in  the 
clay.  However,  due  to  the  presence  of  cracks  in  the  clay,  there  may  be  “enhanced 
diffusion”  of  dissolved  phase  CAH  in  the  clay  due  to  advective  transport  in  the  cracks. 
Advective  transport  in  the  cracks  coupled  to  diffusive  transport  in  the  clay  matrix  is 
simulated  using  a  dual  domain  model.  Model  results  indicate  that  as  the  cracking 
increases,  there  is  an  increase  in  the  mass  of  contaminant  that  enters  and  is  stored  in  the 
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clay.  This  leads  to  higher  concentrations  of  contaminant  that  “back  diffuses”  out  of  the 
clay,  resulting  in  higher  downgradient  plume  concentrations  (known  as  “tailing”)  long 
after  the  contaminant  source  has  been  removed. 
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MODELING  ENHANCED  STORAGE  OF  GROUNDWATER  CONTAMINANTS 
DUE  TO  THE  PRESENCE  OF  CRACKS  IN  LOW  PERMEABILITY  ZONES 
UNDERLYING  CONTAMINANT  SOURCE  AREAS 


1.0  Introduction 


1.1  Background 

Throughout  the  Air  Force  and  the  Department  of  Defense  (DoD),  chlorinated 
solvents  such  as  trichloroethylene  (TCE)  and  tetrachloroethylene  (PCE),  which  are 
classified  as  Dense  Non-Aqueous  Phase  Liquids  (DNAPLs),  have  been  used  as  cleaning 
solutions  to  clean  items  as  diverse  as  aircraft  and  motor  vehicle  parts  and  clothing  (PCE 
is  commonly  used  for  dry-cleaning).  Before  environmental  regulations  were  in  place 
banning  the  improper  disposal  of  these  solvents,  they  were  often  disposed  of  in  sanitary 
and  stonn  sewers,  disposal  pits,  and  fire  training  areas.  Disposal  in  this  manner  has  led  to 
subsurface  soil  and  groundwater  contamination.  DNAPLs  are  denser  than,  and  only 
slightly  soluble  in,  water.  Due  to  their  density,  DNAPLs  move  vertically  through  soils 
and  groundwater  until  their  downward  movement  is  impeded  by  low  penneability  layers. 
The  DNAPL  will  then  pool  atop  the  low  permeability  layers.  Figure  1  is  a  conceptual 
illustration  showing  DNAPL  migration  in  the  subsurface.  Particularly  note  the  DNAPL 
that  has  pooled  atop  the  two  low  penneability  layers. 
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Figure  1:  DNAPL  Movement  in  the  Subsurface  (after  Heiderscheidt,  2010) 


Due  to  the  large  amount  of  chlorinated  solvents  that  the  DoD  and  industry  have  used  and 
disposed  of  over  the  years,  both  properly  and  improperly,  chlorinated  solvent 
contamination  is  an  important  issue.  The  U.S.  Environmental  Protection  agency  (EPA) 
estimates  that  60%  of  Superfund  sites  across  the  United  States  are  contaminated  with 
DNAPLs  (U.S.  EPA,  2010). 

Methods  of  remediating  DNAPL  source  areas,  such  as  surfactant/cosolvent 
flushing  or  in  situ  chemical  oxidation/reduction,  involve  introducing  a  solution  into  the 
subsurface  in  order  to  solubilize,  mobilize,  or  chemically  react  with  the  NAPL.  The 
introduced  solution  flows  through  the  high  permeability  zones,  and  typically  does  not 
interact  with  DNAPL  stored  in  the  low  permeability  layers.  Thus,  even  though  DNAPL 
may  have  been  removed  from  the  high  permeability  zones,  the  low  permeability  layers 
serve  to  store  DNAPL.  Thus  a  concentration  gradient  is  established;  high  concentrations 
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in  the  low  permeability  zones,  and  low  concentrations  in  the  high  permeability  zones, 
which  have  been  remediated.  Because  of  the  concentration  gradient,  DNAPL  can  then 
diffuse  back  from  the  low  penneability  layers  into  the  high  permeability  zones.  This  is 
known  as  “back  diffusion.” 

The  storage  of  DNAPLs  in  low  permeability  lenses  and  layers  in  the  subsurface 
creates  long-tenn  contaminant  sources.  As  this  storage  limits  the  ability  of  remediation 
technologies  to  meet  cleanup  goals,  and  impacts  the  longevity  of  groundwater 
contamination  plumes,  it  is  necessary  that  we  develop  a  better  understanding  of  the 
storage  process.  It  is  thought  that  contaminant  movement  into  these  low  permeability 
lenses  and  layers  occurs  primarily  through  the  process  of  transverse  diffusion  (Ball  et  ah, 
1997;  AFCEE,  2007;  Parker  et  ah,  2008).  However,  there  is  also  evidence  that  either 
advective  transport  or  enhanced  diffusive  transport  occurs,  as  reported  contaminant 
concentrations  in  these  low  permeability  zones  can  be  greater  than  that  which  can  be 
accounted  for  by  diffusion  alone  (AFCEE,  2007).  This  enhanced  transport  may  be  due  to 
naturally  occurring  cracks  present  in  the  low  penneability  zone  or  cracking  due  to 
DNAPL  that  has  pooled  atop  the  low  permeability  zone.  The  interaction  between  the 
DNAPL  and  low  permeability  zone  may  affect  the  properties  of  the  clay  comprising  the 
low  permeability  zones,  thereby  promoting  enhanced  transport  of  the  DNAPL  into  these 
zones  (Ayral  et  ah,  2010).  It  has  been  reported  that  the  hydraulic  conductivity  of 
common  clay  minerals  can  be  one  to  five  orders  of  magnitude  times  greater  when 
penneated  by  organic  liquids  than  when  penneated  by  water,  depending  on  the  organic 
liquid,  and  the  fraction  and  type  of  clay  (Brown  and  Thomas,  1987).  This  increase  in 
hydraulic  conductivity  has  been  ascribed  to  interlayer  compression  (Li  et  al.,  1996). 
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Essentially,  contact  with  organic  liquids  causes  the  clay  structure  to  shrink,  leading  to  the 
formation  of  cracks  and  microfractures  and  a  concomitant  increase  in  hydraulic 
conductivity  (Brown  and  Thomas,  1987).  Current  modeling  approaches  looking  at  the 
transport  of  contaminants  into  and  out  of  clay  layers  and  lenses,  such  as  those  utilized  by 
Ball  et  al.  (1997),  Parker  et  al.  (2008),  and  Sale  et  al.  (2008)  treat  the  process  as  diffusion 
in  a  dilute  solution,  modified  by  a  tortuosity  factor  to  account  for  the  reduced  area  and  the 
increased  transport  path  length  in  a  porous  medium.  This  approach  does  not  account  for 
the  cracks  in  the  clay,  which  may  be  pre-existing  or  the  result  of  interactions  between 
DNAPL  pools  and  low  permeability  layers  and  lenses.  As  discussed  above,  the  impact  of 
these  interactions  may  be  significant  (e.g.,  Brown  and  Thomas,  1987).  Not  accounting 
for  the  altered  clay  structure,  and  the  ensuing  enhanced  transport  into  these  low 
permeability  layers  that  might  result,  could  lead  to  significantly  erroneous  predictions  of 
storage  in  the  low  permeability  zones  and  the  rates  at  which  the  DNAPLs  enter  into  and 
are  released  from  these  zones.  If  the  analysis  of  site  remediation  options  are  then  based 
on  the  faulty  predictions  of  storage  and  transport,  inappropriate  strategies  may  be 
recommended  for  managing  the  site. 

Researchers  at  the  Air  Force  Institute  of  Technology  (AFIT),  in  collaboration  with 
researchers  at  the  University  of  Michigan,  have  proposed  to  model  the  enhanced  transport 
and  storage  of  groundwater  contaminants  in  low  penneability  layers  that  underlie  a 
contaminant  source  area.  The  researchers  at  the  University  of  Michigan  will  focus  on 
experimentally  studying  how  low  permeability  soil  materials  are  modified  in  the  presence 
of  both  pure  DNAPFs  and  DNAPFs  with  surfactants.  It  was  proposed  that  these 
experimental  results  would  then  be  used  by  AFIT  researchers  to  model  the  effect  of 
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DNAPLs  on  low  penneability  materials,  and  ultimately  apply  the  model  to  simulate  (1) 
the  transport  of  DNAPL  compounds  into  and  out  of  the  low  permeability  zones,  and  (2) 
the  evolution  of  dissolved  phase  contaminant  plumes  from  the  low  permeability  source. 

Throughout  this  chapter  and  in  subsequent  chapters  any  reference  to  DNAPL  will 
be  for  the  separate,  immiscible  phase.  The  dissolved  contaminant  phase  will  be  defined 
and  referred  to  as  dissolved  CAH. 

1.2  Research  Objective 

The  primary  objective  of  this  research  is  to  develop  a  model  that  demonstrates 
enhanced  transport  of  groundwater  contaminants  into  and  out  of  a  low  permeability  layer 
from  a  DNAPL  source  atop  the  layer  and  simulate  how  this  enhanced  transport  might 
affect  the  evolution  and  longevity  of  the  dissolved  phase  plume.  The  enhanced  transport 
is  assumed  to  be  the  result  of  cracks  in  the  low  permeability  layer. 

1.3  Specific  Research  Questions 

1 .  What  processes  occur  when  DNAPLs  and  low  permeability  layers  interact? 

2.  What  is  the  effect  of  cracking  on  the  transport  of  DNAPLS  into  and  out  of  low 
permeability  layers? 

3.  What  mathematical  models  are  available  that  can  be  used  to  model  diffusion  and 
enhanced  diffusion  into  low  permeability  layers? 

4.  What  mathematical  models  can  be  used  to  simulate  other  processes  (e.g.,  advection, 
dispersion,  diffusion)  that  have  been  found  to  contribute  to  enhanced  transport  of  the 
DNAPL  into  the  low  permeability  layer? 
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5.  What  is  the  effect  of  enhanced  transport  into  low  permeability  layers  on  dissolved 
plume  longevity  and  evolution? 

1.4  Research  Approach 

The  initial  phase  of  the  study  involves  a  literature  review  to:  1)  determine  the 
appropriate  set  of  mathematical  equations  that  can  be  used  to  model  diffusion  into  low 
penneability  zones;  2)  learn  what  processes  occur  when  DNAPLs  and  low  permeability 
layers  interact;  3)  determine  the  appropriate  set  of  mathematical  equations  that  can  be 
used  to  model  enhanced  diffusion  into  low  permeability  zones. 

The  second  phase  of  the  study  involves  coupling  diffusion  and  enhanced  diffusion 
submodels  with  “standard”  dissolved  contaminant  transport  models  (advection-dispersion 
equations)  and  use  the  coupled  model  to  simulate  the  effect  of  enhanced  diffusion  on 
plume  longevity  and  growth. 

1.5  Scope  and  Limitation  of  Research 

Parameter  values  used  in  the  enhanced  diffusion  model  are  “guesstimated”  based 
on  the  literature.  No  “real  world”  remediation  site  data  are  available  to  validate  the 
model.  Also,  the  model  developed  in  this  study  is  a  simplified  model  that  assumes  (1)  a 
sorbing,  but  otherwise  conservative  contaminant  (no  degradation),  (2)  the  subsurface 
material  properties  are  homogeneous  in  space  and  time;  (3)  steady  state  flow.  Also  note 
that  the  model  simulates  dissolved  contaminant  transport  only;  the  DNAPL  source  is 
assumed  to  be  immobile.  These  assumptions  are  made  so  that  the  effect  of  enhanced 
diffusion  on  plume  evolution  can  be  focused  on.  Finally  the  concentration  versus  time 
breakthrough  data,  which  are  used  to  study  plume  longevity,  are  simulated  at  a  single 
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observation  point  downgradient  from  the  DNAPL  source.  Different  concentration  versus 
time  data  would  be  simulated  if  the  observation  point  were  moved  closer  or  further  from 
the  DNAPL  source.  Nevertheless,  it  is  assumed  that  the  overall  plume  behavior  can  be 
understood,  at  least  in  a  qualitative  sense,  by  observing  the  breakthrough  behavior  at  a 
representative  downgradient  observation  point 

1.6  Definition  of  Terms 

Advection  -  fluid  transport  due  to  fluid  motion 

Basal  spacing  -  spacing  between  adjacent  layers  of  a  crystalline  structure 

Crack  -  an  opening  in  a  material  caused  by  an  applied  stress  that  allows  fluid  to  freely 

enter  the  material 

Dense  Non-aqueous  Phase  Liquid  (DNAPL)  -  a  liquid  that  is  denser  than  water  and  also 
immiscible  in  water. 

Dielectric  Constant  (E)  -  ratio  of  the  permittivity  of  a  material  to  the  permittivity  of  a 
vacuum.  Low  E  materials  have  a  low  ability  to  polarize  and  hold  a  charge  and  thus  have 
low  permittivity. 

Diffusion  -  movement  of  a  dissolved  solute  caused  by  a  concentration  gradient 
Dispersion  -  spreading  of  mass  beyond  a  region  it  normally  would  occupy  due  to 
advection  alone 

Electrophoteric  potential  -  movement  of  particles  under  an  applied  electric  field. 
Enhanced  Diffusion  -  diffusion  from  cracks  in  a  material  that  allows  increased  mass  to 
enter  the  material  due  to  the  presence  of  the  cracks 
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Hydraulic  conductivity  -  the  factor  of  proportionality  between  specific  discharge  and 
hydraulic  gradient 

Permeability-  the  ease  with  which  fluid  can  move  through  porous  material 
Pennittivity  -  measure  of  a  materials  ability  to  be  polarized  by  an  electric  field. 

Zeta  potential  -  a  calculated  value  that  indicates  the  degree  of  repulsion  between  two 
adjacent,  similarly  charged  particles. 
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2.0  Literature  Review 


2.1  Overview 

Groundwater  contaminants,  such  as  chlorinated  aliphatic  hydrocarbons  (CAHs), 
stored  in  subsurface  low  permeability  layers  pose  long  term  challenges  for  remediation. 

It  is  estimated  that  over  3,000  Department  of  Defense  (DOD)  sites  are  contaminated  with 
various  CAHs  which  were  used  as  solvents  (SERDP,  2008).  Dissolved-phase  plumes 
tend  to  be  large  and  can  extend  for  miles  downgradient  from  the  original  contamination 
source.  Contamination  plumes  that  contain  CAHs  are  difficult  to  remediate  due  to  the 
relatively  high  solubility  and  limited  retardation  of  CAH  transport  due  to  sorption 
(SERDP,  2008).  The  ability  of  low  penneability  layers  to  store  CAHs  further  hinders  the 
CAH  contamination  remediation  effort.  In  a  situation  where  high  penneability  layers, 
which  could  serve  as  drinking  water  sources,  are  contaminated  with  CAHs,  and  the 
subsurface  also  consists  of  low  permeability  layers,  these  low  permeability  layers  can  act 
as  CAH  contamination  sources  after  the  high  permeability  zones  are  remediated.  This 
process  of  aquifer  recontamination  is  caused  by  the  back  diffusion  of  the  CAH  from  the 
low  permeability  layer  to  the  high  permeability  layer.  This  back  diffusion  is  caused  by 
the  concentration  gradient  that  is  set  up  between  the  clean  high  permeability  layer  and  the 
contaminated  low  permeability  layer. 

A  better  understanding  of  the  impact  that  chlorinated  solvents  stored  in  low 
penneability  layers  has  is  needed  to  improve  the  ability  to  reach  groundwater  cleanup 
goals.  The  CAHs  that  are  the  focus  of  this  research  are  dense  nonaqueous  phase  liquids 
(DNAPLs)  like  trichloroethylene  (TCE). 
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2.2  Transport  of  Dissolved  Compounds  in  Low  Permeability  Materials 

The  primary  transport  mechanism  for  dissolved  contaminants  in  low  penneability 
layers,  where  advection  is  negligible,  is  simple  Fickian  diffusion.  Johnson  et  al.  (1989) 
used  field  data  collected  from  a  waste  disposal  facility  in  Ontario  to  show  that  transport 
by  simple  Fickian  diffusion  is  more  important  than  advective  transport  when  assessing 
the  transport  of  contaminants  into  low  permeability  layers.  Core  samples  were  collected 
from  an  unweathered  clay  deposit  beneath  a  five  year  old  hazardous  waste  landfill  and 
analyzed  for  chlorides  and  volatile  organic  compounds.  A  plot  of  chloride  concentration 
versus  depth  below  the  landfill/clay  interface  showed  that  the  maximum  depth  that  the 
chloride  was  transported  was  83  cm.  The  plot  showed  the  chloride  concentration 
decreasing  with  distance  and  suggested  chloride  transport  was  due  to  “classic  diffusion” 
(Johnson  et  al.,  1989).  Using  mathematical  modeling  the  authors  were  able  to  confirm 
that  the  dominant  chloride  transport  mechanism  was  diffusion  in  accordance  with 
Equation  2.1,  Fick’s  second  law, 


dc  _  d2C 

—  —  Upff  ~  7 


dt 


'eff  dy2 


(Eq  2.1) 


where  C  is  concentration  [ML'3],  Deff  [L2T_1]  is  the  effective  diffusion  coefficient,  and  y 
[L]  is  the  vertical  coordinate,  with  the  positive  direction  downwards.  The  effective 
diffusion  coefficient  (Deff)  is  related  to  the  free-solution  diffusion  coefficient  (D0)  by  the 
following  equation: 


Deff  =  D0x  (Eq  2.2) 

where  x  is  a  tortuosity  factor  (0<x<l)  that  accounts  for  the  tortuous  nature  of  the  clay 
material  (Clark,  2009). 
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The  organic  compounds  were  only  detected  down  to  a  depth  of  15cm.  The 
concentration  versus  depth  data  for  the  organic  compounds  also  could  be  simulated 
assuming  Fickian  diffusion.  Sorption  of  the  organic  compounds  to  the  clay  particles  was 
assumed  to  play  an  important  role  in  retarding  the  movement  of  the  organic  compounds 
through  the  low  permeability  layer  (Johnson  et  ah,  1989). 

The  effect  of  sorption  of  the  organic  compounds  to  the  clay  particles  can  be 
incorporated  into  the  diffusion  equation  by  modifying  the  effective  diffusion  coefficient 

Drff  =  (Eq  2  3> 

l'imm 

where  Rimm  [-]  is  a  retardation  factor  that  accounts  for  linear,  equilibrium  sorption  of  the 
organic  compound  to  the  clay  immobile  phase: 

Rimm  =  l+£Yi  (Eq  2.4) 

where  pb  [ML'  ]  is  the  bulk  density  of  clay,  0  [-]is  the  porosity  of  clay,  and  Kd  [M'  L  ] 
is  the  sorption  partition  coefficient  of  dissolved  organic  compound  (Clark,  2009). 

Parker  et  al.  (2004)  conducted  a  field  study  at  an  industrial  site  in  Connecticut  to 
detennine  the  integrity  of  a  20-m-thick  clayey  silt  aquitard  that  had  an  accumulated  pool 
of  free-product  trichloroethylene  (TCE)  dense  non-aqueous  phase  liquid  (DNAPL)  atop 
it.  From  1952  to  1972  TCE  was  used  at  the  facility  for  degreasing  operations.  It  was 
assumed  that  residual  and  free-product  TCE  started  pooling  atop  the  aquitard  sometime  in 
the  1950s.  The  sand  aquifer  below  the  aquitard  was  pumped  in  order  to  provide  water  to 
the  site  from  the  1950s  to  2001.  Parker  et  al.  (2004)  proposed  that  TCE  could  be 
transported  into  the  aquitard  by  one  or  more  of  the  following  processes:  (1)  bulk 
movement  of  the  DNAPL  (free-product  TCE),  (2)  advection  of  dissolved  TCE  due  to  the 
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hydraulic  gradient  established  by  pumping  the  lower  aquifer,  and/or  (3)  diffusion  of 
dissolved  TCE.  The  researchers  sampled  for  dissolved  and/or  separate  phase  TCE  in  the 
aquitard,  in  order  to  ascertain  the  primary  transport  processes  within  the  low  permeability 
clayey  silt. 

Cores  collected  from  the  aquitard  indicated  the  presence  of  TCE  at  a  depth  of  2.5 
to  3  m  below  the  top  of  the  aquitard  (Parker  et  ah,  2004).  These  results  were  compared  to 
the  results  of  one-dimensional  analytical  and  numerical  solutions  of  TCE  migration  in  the 
aquitard  which  showed  downward  migration  of  aqueous  TCE  to  a  depth  of  4  to  5  m  into 
the  aquitard  (Parker  et  al.,  2004).  Cases  assuming  diffusion  only  and  diffusion  and 
advection  were  compared.  Based  on  measured  TCE  concentration  values  and  the  results 
of  the  one-dimensional  model  simulations,  the  conclusion  was  that  diffusion  of  aqueous 
TCE  is  the  dominant  transport  process  in  the  clayey  silt  aquitard  and  that  advection,  due 
to  the  hydraulic  gradient  established  by  pumping  the  underlying  aquifer,  was  negligible. 

2.3  Back  Diffusion 

As  stated  previously,  back  diffusion  can  occur  when  a  concentration  gradient  is 
established  between  an  uncontaminated  high  permeability  zone  and  a  contaminated  low 
penneability  zone.  Parker  et  al.  (2008)  showed  that  back  diffusion  from  an  aquitard  in  a 
sand  aquifer  could  cause  contaminant  persistence  above  MCLs  long  after  the  contaminant 
plume  source  zone  was  isolated  or  removed. 

Wilson  (1997)  studied  how  pump-and-treat  systems  used  to  remediate 
contaminated  aquifers  have  failed  to  restore  aquifers  to  drinking  water  standards. 

Aquifer  heterogeneities  (i.e.,  numerous  high  and  low  permeability  zones  interspersed  in 
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the  aquifer)  were  identified  as  one  possible  explanation  why  pump-and-treat  systems  fail. 
Contaminant  concentrations  at  pumping  wells  used  to  capture  plumes  typically  drop  to 
low  levels,  where  they  level  off.  On  concentration  versus  time  breakthrough  curves 
constructed  from  measurements  of  concentration  in  water  extracted  at  the  pumping  well, 
this  appears  as  a  long-lasting  low  concentration  tail.  When  pumping  is  stopped  for  a  time 
and  then  resumed,  concentration  levels  are  often  found  to  have  risen  (so-called  rebound). 
Tailing  and  rebound  of  contaminant  concentrations  at  extraction  wells  during  pump-and- 
treat  remediation  are  often  attributed  to  slow  back  diffusion  of  stored  contaminant  from 
low  permeability  zones  within  the  aquifer  (Wilson,  1997). 

Liu  and  Ball  (2002)  studied  the  diffusion  of  tetrachloroethene  (PCE)  and  TCE  in  a 
natural  aquitard  below  Dover  Air  Force  Base,  Delaware.  At  this  Dover  Air  Force  Base 
site  bentonite-sealed  steel  sheetpiles  were  installed  to  establish  two  experimental  test 
cells  that  hydraulically  isolated  the  contaminated  aquifer  and  aquitard  from  the 
surrounding  groundwater  flow.  The  sheetpiles  were  installed  to  a  depth  of  approximately 
2.0  m  into  the  aquitard.  The  aquifer  and  aquitard  had  been  exposed  to  water 
contaminated  with  dissolved  (but  not  DNAPL)  PCE  and  TCE  for  more  than  a  decade 
(Mackay  et  ah,  2000).  Three  injection  and  three  extraction  wells  were  installed  in  each  of 
the  two  cells.  One  test  cell  was  subjected  to  continuous  pumping  and  the  other  test  cell 
was  subjected  to  pulsed  pumping.  Over  a  35  month  period  the  test  cells  were  subjected 
to  no-flow  conditions  for  11.1  months  followed  by  8.1  months  of  pumping  followed  by 
15.8  months  of  no-flow  conditions.  Throughout  the  35  month  period,  core  samples  were 
collected  from  the  test  cells  and  analyzed  for  PCE  and  TCE  content.  The  pumping 
conducted  within  the  two  test  cells  and  the  no-flow  conditions  following  the  pumping 
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established  the  necessary  conditions  to  promote  back  diffusion  of  contaminants  from  the 
aquitard  into  the  aquifer  (Liu  and  Ball,  2002). 

Sorption  and  diffusion  data  were  obtained  in  the  laboratory.  These  data  were  used 
in  a  model  that  assumed  sorption-retarded  diffusion  of  TCE  and  PCE  within  the  aquitard. 
Model  simulations  were  run  with  these  independently  detennined  parameters.  The  model 
predicted  concentrations  were  compared  with  coring  results  and  close  agreement  was 
found.  The  authors  concluded  that  diffusion,  not  advection  or  transformation,  was  the 
dominant  transport  process  occurring  in  the  aquitard  (Liu  and  Bah,  2002). 

Chapman  and  Parker  (2005)  investigated  how  back  diffusion  from  a  clayey  silt 
aquitard  caused  the  persistence  of  a  TCE  plume  330  m  downgradient  from  the  TCE 
source,  even  though  the  source  itself  was  isolated  in  1994  by  the  installation  of  a  sheet 
piling  enclosure.  This  investigation  took  place  at  the  same  industrial  site  in  Connecticut 
that  Parker  et  al.  (2004)  reported  on  in  their  investigation  of  TCE  diffusion  in  the  clay 
aquitard.  Data  from  aquitard  cores  helped  detennine  the  contaminant  distribution  in  the 
aquitard  below  the  plume  and  groundwater  data  collected  from  a  site-wide  network  of 
conventional  monitoring  wells  helped  determine  the  concentration  of  TCE  in  the 
groundwater  downgradient  from  the  DNAPL  source. 

Seven  to  ten  months  after  the  DNAPL  source  was  isolated,  data  from  two  of  the 
monitoring  wells  330  m  down  gradient  from  the  isolated  DNAPL  source  indicated  a 
substantial  decrease  in  TCE  concentrations.  The  TCE  concentration  decreased  from 
between  5000  to  30,000  ug/L  to  between  200  to  3000  pg/L  2.5  years  after  installation  of 
the  sheet  pile  enclosure. 
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The  time  between  the  sheet  piling  installation  to  isolate  the  source  zone  and  the 
TCE  concentration  decline  at  the  downgradient  monitoring  well  was  in  close  agreement 
to  the  travel  times  determined  using  the  velocity  calculated  by  Darcy’s  Law  and 
measured  by  borehole  dilution.  This  agreement  confirms  that  the  TCE  concentration 
decline  was  caused  by  the  isolation  of  the  DNAPL  source  zone  (Chapman  and  Parker, 
2005).  Due  to  the  low  organic  carbon  content  of  the  aquifer  sand,  TCE  retardation  was 
assumed  negligible. 

A  decade  after  the  DNAPL  source  isolation,  a  “tail”  of  persistent,  low  TCE 
concentrations  were  measured  in  the  two  monitoring  wells  downgradient  from  the  sheet 
pile  enclosure.  The  conclusion  was  that  sources  of  TCE  outside  of  the  enclosure  existed. 
Three  hypotheses  about  the  source  of  the  TCE  were  proposed.  The  first  hypothesis  was 
that  back  diffusion  of  TCE  from  the  aquitard  below  the  plume  was  occurring.  The 
second  hypothesis  was  that  a  small  area  of  DNAPL  contamination  which  had  been  found 
outside  the  enclosure  was  the  source  of  the  TCE.  The  last  hypothesis  was  that  there  was 
some  other  unknown/unidentified  DNAPL  source  zone  outside  the  enclosure.  Numerical 
modeling  was  used  to  show  that  back  diffusion  fully  accounted  for  the  plume  tailing 
(Chapman  and  Parker,  2005). 

Similar  to  the  study  presented  in  Parker  et  al.  (2004)  of  the  industrial  site  in 
Connecticut,  Parker  et  al.  (2008)  is  an  investigation  of  an  industrial  site  in  Florida  that 
had  an  aquifer  contaminated  with  TCE.  This  site  was  near  a  metal  fabricating  and 
precision  cleaning  facility  that  used  TCE  as  a  solvent.  The  site’s  TCE  contamination  was 
estimated  to  have  occurred  from  the  mid-  to  late  1960s  to  1977.  The  aquifer  below  the 
site  was  described  as  heterogeneous  with  beach  sand  deposits  and  thin  laterally 
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discontinuous  silt  and  clay  lenses.  A  5  to  20  cm  thick  clay  layer  runs  continuously  below 
the  site  at  a  depth  of  8-10  m  below  ground  surface  (bgs).  This  clay  layer  extends  beyond 
the  property  boundary  in  the  direction  of  groundwater  flow  and  separates  the  aquifer  into 
an  upper  and  lower  zone.  The  TCE  contamination  plume  was  caused  by  TCE  DNAPL 
lenses  resting  above  and  below  the  clay  layers  (Parker  et  al.,  2008) 

The  TCE  DNAPL  source  zone  was  isolated  in  2002  using  a  groundwater 
extraction,  treatment,  and  re-injection  system.  The  extraction  and  re-injection  system 
was  intended  to  create  a  clean  water  front  to  flush  the  downgradient  plume.  Though  the 
containment  system  was  in  operation  since  2002,  detectable  levels  of  TCE  and 
degradation  products  still  were  measured  downgradient  of  the-injection  wells.  The 
researchers  proposed  three  hypotheses  on  why  the  plume  persisted  downgradient  of  the 
re-injection  wells.  The  first  hypothesis  was  incomplete  source-zone  capture.  The  second 
hypothesis  was  the  occurrence  of  DNAPL  downgradient  of  the  re-injection  wells.  The 
third  hypothesis  was  that  there  was  back  diffusion  from  the  thin  clay  layers  into  the  sand 
aquifer.  The  first  two  hypotheses  were  eliminated,  leaving  the  explanation  that  back 
diffusion  from  the  clay  layers  caused  the  dissolved  TCE  plume  to  persist.  To  support  the 
hypothesis  that  back  diffusion  was  the  cause  of  the  plume  persistence,  the  researchers 
used  VOC  concentration  data  collected  from  within  and  near  the  clay  layers  and 
numerical  model  simulations  (Parker  et  al.,  2008).  Their  analysis  is  described  in  detail, 
below. 

Before  source  isolation,  groundwater  sampling  was  conducted  to  detennine  the 
contaminant  distribution  downgradient  of  the  source  zone.  The  highest  contaminant 
concentrations  occurred  close  to  the  clay  layers.  After  the  source  isolation  groundwater 
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samples  were  once  again  collected.  The  results  of  those  samples  showed  high 
contaminant  concentrations  in  the  vicinity  of  the  clay  layers  even  after  contaminant 
concentrations  would  have  been  expected  to  decline  due  to  isolation  of  the  source.  Core 
samples  collected  indicated  VOC  mass  storage  in  the  clay  layers. 

Using  numerical  modeling,  simulations  were  run  to  show  if  back  diffusion  is 
capable  of  producing  the  observed  field  results  (Parker  et  al.,  2008).  The  numerical 
modeling  was  conducted  using  HydroGeoSphere.  Simulations  were  run  with  a  single 
suspended  clay  layer  in  the  sand  aquifer  and  with  several  thin,  discontinuous  clay  layers 
along  with  an  underlying  aquitard.  For  the  case  of  the  single  suspended  clay  layer  the 
model  was  2-D  with  uniform  flow  occurring  in  two  aquifers  separated  by  a  0.2  m  thick 
clay  layer.  The  contaminant  source  was  a  thin  TCE  DNAPL  layer  along  the  top  of  the 
clay  near  the  upgradient  end  of  the  model  domain.  The  concentration  of  the  TCE  in  the 
water  adjacent  to  contaminant  source  was  the  TCE  solubility  limit  and  the  source  was 
assumed  to  have  accumulated  over  a  35-year  period.  After  the  35-year  period,  the  source 
concentration  was  set  to  zero  which  allowed  for  clean  water  to  flush  through  the  system. 
This  simulation  was  intended  to  isolate  the  back  diffusion  effects  from  a  single  clay  layer. 
The  results  of  the  simulation  showed  that  50  years  after  the  source  is  removed  the  plume 
still  lingered  (Parker  et  al.,  2008). 

A  simulation  was  also  run  where  the  aquifer  contained  several  suspended,  thin, 
discontinuous  clay  layers  along  with  an  underlying  clay  aquitard.  The  thickness  of  these 
clay  layers  was  0.4  m  to  0.6  m.  For  a  30-year  period,  a  constant  concentration  of  TCE  at 
its  solubility  limit  was  placed  on  top  of  the  clay  layers.  At  the  end  of  the  30-year  period, 
the  source  was  removed  and  the  aquifer  was  flushed  with  clean  water  for  200  years.  The 
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plume  still  persisted  50  years  after  the  source  was  removed.  After  200  years, 
concentrations  at  some  locations  in  the  model  dropped  to  below  0.01  mg/L  except  in  the 
vicinity  of  “one  thicker  clay  layer”  (Parker  et  al.,  2008).  In  the  vicinity  of  the  clay 
aquitard  at  the  bottom  of  the  aquifer,  higher  contaminant  concentrations  still  persisted.  It 
was  noted  that  diffusion  out  of  the  clay  aquitard  was  what  contributed  to  these  higher 
concentrations.  Parker  et  al.  (2008)  concluded  that  the  stored  mass  in  the  clay  aquitard 
was  released  more  slowly  as  compared  to  the  release  from  the  suspended  clay  layers 
because  only  the  top  surface  of  the  clay  aquitard  is  flushed  and  the  contaminant  diffused 
deeper  into  the  clay  aquitard  compared  to  the  suspended  clay  layers. 

Sale  et  al.  (2008)  studied  the  effect  that  back  diffusion  from  low  permeability 
layers  had  on  downgradient  water  quality.  A  two-layer  system  consisting  of  a 
transmissive  sand  layer  and  a  low  permeability  silt  layer  was  studied.  In  the  transmissive 
layer,  a  high  concentration  source  analogous  to  a  thin  DNAPL  pool  was  introduced  above 
the  low  penneability  layer.  The  source’s  dissolved  constituents  were  transported  in  the 
transmissive  zone  by  advection  and  by  transverse  diffusion  in  the  low  penneability  layer. 
After  the  source  was  removed,  the  concentration  gradient  between  the  low  permeability 
layer  and  the  transmissive  layer  was  reversed  causing  back  diffusion  of  the  contaminants 
from  the  low  permeability  layer  into  the  transmissive  layer.  This  back  diffusion  can  lead 
to  a  sustained  contaminant  concentration  in  the  transmissive  layer  after  the  removal  of  the 
contaminant  source.  The  authors  noted  that  back  diffusion  occurs  over  an  even  longer 
time  scale  than  the  initial  diffusion  into  the  clay  layer,  with  important  implications 
regarding  site  remediation  (Sale  et  al.,  2008). 
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2.4  DNAPL  -  Low  Permeability  Material  Interaction 

DNAPL  in  contact  with  low  penneability  layer  material  may  alter  and  modify  the 
material’s  physical  properties,  hence  allowing  enhanced  dissolved  phase  CAH  transport 
into  the  low  permeability  layers.  Brown  and  Thomas  (1987)  studied  the  mechanism  by 
which  low  dielectric  constant  (E)  organic  chemicals  increase  the  hydraulic  conductivity 
of  compacted  clay  materials.  Hydraulic  conductivity  measurements,  basal  spacing, 
electrophoretic  mobility,  zeta  potential,  flocculation,  and  volume  change  were  measured 
for  clay  minerals,  kaolinite,  illite,  and  smectite  (Brown  and  Thomas,  1987).  Since  the 
focus  of  this  current  research  is  on  the  interaction  between  DNAPLs  and  low 
permeability  clay  materials,  the  results  of  the  measurements  of  hydraulic  conductivity, 
basal  spacing  and  volume  change  that  Brown  and  Thomas  (1987)  discuss  will  be 
presented  in  this  section. 

The  hydraulic  conductivity  measurements  were  made  by  compacting  mixtures  of 
kaolinite  and  water,  illite  and  water  and  smectite  and  water  in  10-cm  diameter  fixed  wall 
penneameters.  Replicates  of  each  mixture  were  exposed  to  solutions  of  “0,  60,  80,  and 
100%  (by  volume)  ethanol  in  distilled  deionized  water;  0,  60,  80,  and  100%  (by  volume) 
acetone  in  distilled  deionized  water;  and  0.10,  0.20,  and  0.30  M  NaCl”  (Brown  and 
Thomas,  1987).  Results  showed  that  for  E  values  <40,  the  hydraulic  conductivity 
increased  for  kaolinite  penneated  by  acetone.  When  permeated  with  ethanol  the  kaolinite 
hydraulic  conductivity  increased  when  the  solution’s  E  value  was  <45.  Both  the  acetone 
and  ethanol  solutions  increased  the  illite  hydraulic  conductivity  when  the  solution’s  E 
value  was  less  than  45.  A  decreasing  E  value  corresponds  to  an  increasing  organic 
fraction  and  a  decreasing  water  fraction  of  a  solution  (Brown  and  Thomas,  1987). 
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Measurements  of  smectite  clay  basal  spacing  versus  mixtures  of  acetone  and 
ethanol  with  varying  E  values  were  made  using  x-ray  diffraction.  Basal  spacing  is  a 
measurement  of  space  between  adjacent  layers  in  a  crystalline  structure.  In  general,  as 
basal  spacing  increases,  hydraulic  conductivity  decreases.  According  to  Brown  and 
Thomas  (1987)  measurements  of  basal  spacing  were  only  made  for  smectite  clay  because 
only  smectite  clay  undergoes  changes  in  basal  spacing  when  permeated  by  concentrated 
organic  liquids.  The  basal  spacing  for  smectite  equilibrated  with  water  was  1.8  nm. 

Dilute  concentrations  of  acetone  (2-5%)  with  E  values  of  77  and  76,  respectively,  caused 
the  basal  d  spacing  of  the  smectite  to  increase  to  2.0  nm.  An  acetone  in  water  E  value  of 
49  caused  the  basal  d  spacing  to  decrease  to  1.8  nm.  The  basal  spacing  decreased  to  1.4 
nm  when  the  solutions  E  value  was  decreased  to  43  (Brown  and  Thomas,  1987). 

In  dilute  ethanol  solutions  the  smectite  basal  spacing  increased  to  2.3  nm.  For  a 
solution  E  value  of  57  the  spacing  decreased  to  2.0  nm.  For  solution  E  values  below  35 
the  spacing  decreased  to  1.6  nm  (Brown  and  Thomas,  1987). 

The  results  of  the  measurements  showed  that  dilute  concentrations  of  organic 
liquids  increase  the  basal  spacing  of  the  clay  material  which  in  turn  decreases  the 
hydraulic  conductivity  of  the  clay  material  (Brown  and  Thomas,  1987). 

Finally  Brown  and  Thomas  (1987)  measured  the  volume  changes  of  the  clay 
materials  when  exposed  to  acetone,  xylene,  and  water.  Their  results  showed  that  xylene 
(E=2.4)  caused  the  materials  to  swell  the  least  and  water  (E=78)  caused  the  materials  to 
swell  the  most.  The  smectite  material  swelled  the  most,  which  the  authors  attributed  to 
the  increased  basal  spacing.  Changes  in  the  particle  spacing  are  what  caused  the  kaolinite 
and  illite  to  swell  (Brown  and  Thomas,  1987). 
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Brown  and  Thomas  (1987)  came  to  the  conclusion  that  as  organic  liquids  displace 
water  that  is  in  equilibrium  with  clay  soil,  the  soil  will  shrink,  causing  cracks  to  form 
which  can  act  as  channels  for  liquids  to  flow.  This  results  in  an  increase  in  hydraulic 
conductivity. 

Clay-DNAPL  interaction  may  not  be  the  only  cause  of  crack  formation  in  clay. 
Fookes  and  Denness  (1969)  theorized  the  natural  existence  of  cracks  in  clay.  The 
research  performed  by  Fookes  and  Denness  (1969)  was  conducted  to  detennine  the 
fissure  patterns  in  soft  rocks  so  that  the  engineering  behavior  of  fissured  sediments  could 
be  analyzed.  The  fissure  patterns  of  Chalk,  Gault,  and  Weald  Clays  were  examined  and  it 
was  concluded  that  differences  between  fissure  patterns  could  be  explained  by 
considering  the  geology  of  where  the  materials  came  from. 

This  current  research  is  intended  to  develop  a  better  understanding  of  the  impact 
that  pooled  DNAPLs  have  on  low  penneability  clay  layers  and  understand  how  the 
transport  and  storage  of  dissolved  CAHs  are  affected  by  this  impact,  or  alternatively, 
affected  by  the  natural  presence  of  cracks  in  the  clay  matrix.  Brown  and  Thomas  (1987) 
showed  that  organic  liquids  can  increase  the  hydraulic  conductivity  of  clay  materials  by 
shrinking  the  clay,  which  allows  cracks  to  develop,  thereby  allowing  organic  liquids  to 
move  more  freely  into  the  clay  material.  Fookes  and  Denness  (1969)  showed  that 
cracking  is  commonly  found  in  clays.  The  next  section  is  intended  to  establish  the 
framework  for  how  modeling  can  be  used  to  simulate  the  transport  of  dissolved  CAHs 
within  the  cracks  of  low  penneability  clay  layers. 
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2.5  Modeling  Contaminant  Transport  in  Cracks 

Grisak  and  Pickens  (1980)  used  numerical  modeling  to  describe  solute  transport 
through  fractured  media  that  accounted  for  the  effect  of  matrix  diffusion  on  transport. 
Their  conceptual  model  included  mobile  water  in  highly  penneable  fractures  and 
immobile  water  in  a  low  permeability  matrix.  Their  model  accounted  for  the  following: 
1)  diffusion  of  solute  from  the  fractures  into  the  low  permeability  matrix;  2)  advection 
and  dispersion  due  to  water  flow  in  the  fractures;  and  3)  linear,  equilibrium  sorption  in 
the  matrix. 

The  dominant  transport  mechanism  in  a  fracture  is  advection-dispersion,  while  in 
the  unfractured  matrix,  the  dominant  transport  mechanism  is  Fickian  diffusion.  Matrix 
diffusion  is  considered  a  much  slower  transport  process  than  the  advection-dispersion 
transport  process.  Figure  2  shows  solute  transport  in  a  single  fracture  surrounded  by  an 
unfractured  matrix.  The  input  concentration  of  the  solute  is  C0  and  the  fracture  aperture 
width  is  2b.  The  solute  is  being  transported  down  the  fracture  by  advection/dispersion 
and  transported  into  the  surrounding  matrix  by  diffusion.  The  solute  profile  in  the 
fracture  and  the  matrix  is  depicted  by  the  hatched  area  of  the  figure.  There  is  a  zero 
concentration  gradient  (5C/3y=0)  at  the  center  of  each  matrix  block,  which  assumes  that 
there  are  other  fractures  at  a  distance  2B  from  the  center  of  the  fracture  being  simulated. 
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FRACTURE 

•  X 

-z 

Figure  2:  Solute  Transport  in  a  Fracture  of  aperture  2b  (after  Grisak  and  Pickens,  1980) 

Grisak  and  Pickens  (1980)  conducted  sensitivity  analyses  to  determine  how 
various  parameter  values  affected  solute  transport  in  a  fracture  system.  To  study  the 
effect  of  diffusion  within  the  matrix,  Grisak  and  Pickens  (1980)  applied  a  constant 
concentration  source  (concentration  =  Co)  at  z  =  0,  and  plotted  relative  concentration, 
C/C0,  versus  time  breakthrough  curves  at  z  =  -0.76  m  for  matrix  diffusion  coefficients 
ranging  from  0.0  cm  /s  to  10'  cnr/s.  The  breakthrough  curves  indicated  that  as  the 
diffusion  coefficient  got  larger  the  effect  of  matrix  diffusion  became  more  pronounced 
and  the  shallower  the  breakthrough  curves  became.  The  authors  theorized  that  if  the 
concentration  source  was  cut  off,  then  a  concentration  gradient  between  the  fracture  and 
the  matrix  would  be  established  and  the  solute  stored  in  the  matrix  would  be  transported 
by  back  diffusion  into  the  fracture  (Grisak  and  Pickens,  1980). 

Next  the  effect  of  aperture  size  was  considered.  Assuming  a  constant  fluid 
velocity  in  the  fracture,  reducing  the  aperture  size  reduces  the  quantity  of  solute 
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transported  in  the  fracture  and  increases  the  relative  amount  of  solute  that  enters  the 
matrix,  since  for  a  given  diffusion  coefficient,  the  mass  flux  of  solute  into  the  matrix  is 
controlled  by  the  concentration  gradient  only.  When  the  aperture  size  is  decreased  the 
fraction  of  solute  transported  in  the  fracture  is  decreased  while  the  fraction  of  solute 
diffused  into  the  matrix  is  increased.  The  results  are  that  a  relatively  larger  mass  of 
solute  is  retained  near  the  source  at  z  =  0  (Grisak  and  Pickens,  1980). 

To  study  the  effect  of  fluid  velocity  in  the  fracture,  the  velocity  was  varied  from  1 
to  60  m/d  with  other  parameters  held  constant.  The  conclusion  is  that  an  increase  in  fluid 
velocity  leads  to  an  increase  in  the  ratio  of  solute-transported  to  solute-stored  in  the 
matrix  (Grisak  and  Pickens,  1980). 

To  study  the  effect  of  dispersion  in  the  fractures,  the  longitudinal  dispersivity  was 
varied.  The  results  show  that  the  larger  the  dispersivity,  the  earlier  the  breakthrough  at  y 
=  0.76  m.  With  small  values  of  dispersivity,  no  significant  breakthrough  occurs  because 
without  spreading,  there  is  a  large  concentration  gradient  maintained  between  solute  in 
the  fractures  and  the  matrix,  and  therefore,  an  increase  in  solute  diffusion  into  the  matrix. 

The  effect  of  the  matrix  porosity  was  evaluated  for  various  fracture  aperture  sizes. 
The  matrix  porosities  were  varied  and  results  show  that  for  small  matrix  porosities  and 
larger  aperture  sizes  the  solute  concentration  more  rapidly  approaches  its  maximum 
breakthrough  concentration.  This  intuitively  makes  sense,  since  small  matrix  porosities 
result  in  smaller  amounts  of  solute  transported  from  the  fracture  into  the  matrix  (Grisak 
and  Pickens,  1980). 

The  distribution  coefficient,  IQ,  describes  linear,  equilibrium  sorption  of  solute 
between  sorbed  and  dissolved  phases.  Grisak  and  Pickens  (1980)  assumed  sorption  in  the 
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matrix  only.  Varying  the  distribution  coefficient  values,  they  found  that  for  larger  Kd 
values  there  were  larger  net  solute  fluxes  into  the  matrix  from  the  fracture  (Grisak  and 
Pickens,  1980).  This  is  because  that  the  larger  the  distribution  coefficient,  the  less  solute 
is  in  the  dissolved  phase  in  the  matrix,  and  the  higher  the  concentration  gradient  between 
the  fractures  and  matrix. 

Kueper  and  McWhorter  (1991)  studied  the  conditions  that  allow  DNAPL  to  enter 
an  initially  water  saturated  fracture.  They  considered  the  height  of  the  DNAPL  pool  that 
must  accumulate  over  the  fracture  before  the  DNAPL  will  enter  the  fracture  and  used 
numerical  modeling  to  demonstrate  DNAPL  transport  through  fractures.  The  modeling 
shows  that  DNAPL  entry  pressure,  which  is  detennined  by  the  height  of  the  DNAPL  pool 
above  the  fracture,  is  inversely  proportional  to  fracture  aperture  and  directly  proportional 
to  the  interfacial  tension  between  the  DNAPL  and  water  in  the  fracture  (Kueper  and 
McWhorter,  1991). 

Numerical  modeling  also  indicated  that  DNAPL  transport  time  through  a 
fractured  aquitard  is  inversely  proportional  to  the  DNAPL  pool  height,  the  fracture 
aperture  size,  and  the  dip  of  the  fracture  measured  from  the  horizontal  plane.  Downward 
water  gradients  across  the  aquitard  increased  the  DNAPL  flux  into  aquitard  fractures  and 
upward  water  gradients  decreased  the  rate  of  DNAPL  movement  downward  (Kueper  and 
McWhorter,  1991). 

Parker  et  al.  (1994)  developed  a  revised  conceptual  model  of  DNAPL  fate  and 
transport  in  fractured  clay  deposits  and  sedimentary  rock.  They  recognized  the  fact  that 
the  DNAPL  in  the  fracture  can  be  transformed  from  the  immiscible  phase  to  the  dissolved 
phase,  which  then  can  diffuse  into  the  porous  matrix  and  remain  in  the  dissolved  phase  or 
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sorb  onto  the  matrix  material.  In  essence,  Parker  et  al.  (1994)  considered  both  the 
dissolved  transport  represented  in  the  Grisak  and  Pickens  (1980)  model,  and  the  separate 
phase  transport  simulated  by  Kueper  and  McWhorter  (1991).  Factors  identified  that 
affect  transport  are  the  physical  and  chemical  properties  of  the  DNAPL  and  geologic 
material  and  the  extent  to  which  the  DNAPL  is  replenished  from  DNAPL  pools  above  the 
fracture(s)  or  other  DNAPL  filled  fractures.  Moderate  to  narrow  aperture  size,  high 
aqueous  solubility,  large  porosity,  and  high  sorption  capacity  of  the  matrix  were  shown  to 
enhance  the  rate  of  the  immiscible  phase  loss  (Parker  et  al.,  1994). 

Other  published  studies  of  contaminant  transport  through  a  fractured  have  come 
up  with  conclusions  similar  to  the  studies  that  have  been  discussed.  Rowe  and  Booker 
(1990)  concluded  that  in  a  situation  where  a  clay  aquitard  contains  fractures  and  is 
located  above  an  aquifer,  water  quality  degradation  in  the  aquifer  could  be  significant  as 
a  result  of  contaminant  transport  down  through  the  aquitard  fractures  and  into  the  aquifer. 
Harrison  et  al.  (1992)  used  modeling  to  show  that  vertical  fractures  that  fully  penetrate  an 
aquitard  and  have  fracture  aperture  sizes  of  at  least  10  pm  can  provide  an  advective 
transport  pathway  for  dissolved  contaminants  that  can  lead  to  considerable  degradation  of 
water  quality  if  an  aquifer  is  present  below  the  aquitard  (Harrison  et  al.,  1992).  Reynolds 
and  Kueper  (2002)  concluded  that  the  fracture  aperture  is  the  most  important  factor 
controlling  DNAPL  migration  through  a  fracture.  By  increasing  the  fracture  aperture 
from  15  pm  to  50  pm  there  was  a  20  fold  increase  in  the  rate  of  DNAPL  migration 
downward  (Reynolds  and  Kueper,  2002).  They  further  concluded  that  DNAPL  migration 
rates  are  not  significantly  influenced  by  varying  the  parameters  governing  matrix 
diffusion. 
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Models  such  as  the  one  presented  by  Grisak  and  Pickens  (1980)  may  be  useful  in 
simulating  dissolved  solute  transport  into  and  out  of  low  penneability  clay  materials  that 
contain  fractures,  referred  to  as  cracks  from  this  point  on.  This  research  assumes  that  the 
interaction  between  the  DNAPL  and  the  low  penneability  clay  causes  the  cracks  to  form 
in  the  clay,  or  alternatively,  that  cracks  pre-existed  in  the  clay.  The  DNAPL  dissolves, 
and  the  dissolved  CAH  is  transported  by  advection  into  the  cracks.  From  the  cracks,  the 
dissolved  CAH  diffuses  into  the  surrounding  clay  matrix.  The  dissolved  phase  is  also 
transported  in  the  aquifer  by  advection/dispersion  downgradient  of  the  DNAPL  source. 
Chapter  3  will  discuss  the  methodology  used  to  model  these  coupled  transport  processes 
in  order  to  quantify  the  impact  that  enhanced  transport  into  and  out  of  low  penneability 
layers  has  on  dissolved  CAH  plume  evolution  and  longevity. 
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3.0  Methodology 


3.1  Overview 

In  this  chapter,  the  approach  taken  to  model  the  enhanced  transport  of  dissolved 
CAH  into  and  out  of  a  low  penneability  clay  layer  due  to  a  DNAPL  pool  sitting  atop  the 
layer  will  be  presented.  The  model  is  demonstrated  using  the  physical  and  chemical 
properties  of  a  common  CAH  of  particular  importance  to  the  DoD,  TCE.  The 
conceptual  model  consists  of  a  high  conductivity  sand  layer  on  top  of  a  low  conductivity 
clay  layer  (see  Figure  3).  Cracks  are  hypothesized  to  exist  in  the  low  conductivity  clay, 
either  due  to  clay-DNAPL  interaction  or  the  natural  existence  of  cracks  in  the  clay 
(Fookes  and  Denness,  1969).  We  may  model  the  cracks  in  the  clay  as  regions  of  mobile 
water,  while  the  surrounding  non-cracked  clay  matrix  may  be  represented  as  an  immobile 
region.  As  the  cracks  form,  the  mobile  porosity  of  the  clay  increases  as  the  immobile 
porosity  decreases.  In  this  study,  model  simulations  are  run  to  simulate  the  effect  of 
cracking  in  the  clay  layer  on  dissolved  contaminant  transport.  We  hypothesize  that 
“enhanced  diffusion”  due  to  the  cracks  in  the  clay,  which  pennits  contaminant  advection 
through  the  cracks  in  addition  to  diffusion  in  the  clay  itself,  may  significantly  affect 
overall  contaminant  transport. 

3.2  General  Model  Description 

The  key  assumption  of  this  research  is  that  low  permeability  layers  contain 
cracks,  either  naturally  occurring  or  the  result  of  DNAPF  pooling  on  top  of  the  layer. 
Upon  the  formation  of  the  cracks,  the  pooled  DNAPF  dissolves,  and  the  dissolved  CAH 
is  then  transported  into  the  cracks  by  advection  and  dispersion.  From  the  cracks,  the 
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CAH  diffuses  into  the  surrounding  low  permeability  layer.  As  stated  above,  this 
transport  by  advection/dispersion  into  the  cracks  and  diffusion  from  the  cracks  into  the 
surrounding  low  penneability  layer  results  in  “enhanced  diffusion.”  The  assumption  is 
that  this  enhanced  diffusion  process  leads  to  increased  storage  of  CAH  as  compared  to 
the  CAH  storage  that  would  result  from  transport  by  only  Fickian  diffusion  within  the 
intact  low  penneability  layer. 

Figure  3  depicts  a  general  conceptual  model  for  this  research.  A  pool  of  DNAPL 
collects  on  the  surface  of  cracked  low  permeability  clay.  From  the  DNAPL  pool  aqueous 
phase  CAH  is  transported  into  the  cracks  by  advection/dispersion  and  from  the  cracks 
mass  is  transferred  to  the  immobile  region  of  the  low  penneability  clay  by  diffusion. 
Aqueous  phase  CAH  is  also  transported  downgradient  in  the  sand  layer  by  advection, 
dispersion,  and  linear,  equilibrium,  sorption. 
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Figure  3:  Research  Conceptual  Model  Depicting  Advective/Dispersive  Transport  in  Crack, 
and  Mass  Transfer  to  Immobile  Region  by  Diffusion  (represented  as  first  order  mass 

transfer  process) 


Figure  4  depicts  a  simplified  conceptual  model  that  we  can  use  to  study  the 
impact  of  enhanced  diffusion  on  contaminant  transport.  In  essence,  the  scenario  depicted 
in  Figure  4  couples  the  conventional  model  of  advective-dispersive-sorptive  transport  in  a 
high  conductivity  sand  layer  with  the  Grisak  and  Pickens  (1980)  model  of 
advective/dispersive  transport  in  fractures  and  diffusive  transport  in  the  surrounding 
unfractured  low  conductivity  matrix.  In  Figure  4,  dissolved  CAH  from  the  DNAPL  pool 
is  assumed  to  move  into  cylindrical  cracks  of  radius  b  by  advection  and  then  diffuse  from 
the  cracks  into  the  low  permeability  layer.  Dissolved  CAH  also  diffuses  directly  into  the 
low  permeability  layer  from  the  DNAPL  pool.  In  the  high  conductivity  sand,  dissolved 
CAH  is  transported  downgradient  in  the  x-direction  by  advection  and  dispersion.  This 
model  assumes  that  a  linear,  reversible,  equilibrium  relationship  describes  sorption  of 
dissolved  CAH  onto  both  the  sand  and  clay  (both  at  the  cracks  and  within  the  clay 
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matrix),  but  the  extent  of  sorption  on  sand  is  much  less  than  the  extent  of  sorption  on 
clay. 


Figure  4:  Low  Permeability  Layer  with  Cracks  Formed  by  Pooled  DNAPL  Below  a  High 

Permeability  Sand  Aquifer 


3.3  Model  Equations 

Equation  3.1  (Zheng  and  Wang,  1999)  is  the  governing  equation  that  describes 
contaminant  fate  and  transport  in  the  sand  aquifer  layer  above  the  low  permeability  clay. 

R0^  =  A(W^_  A  (0vsC)  (Eq  3.1) 

St  Sxj  V  1J  Sxj/  Sxj  v  1  J  v  U  7 

where  R  [-]  is  the  retardation  factor  for  the  sand  aquifer  layer,  9  [-]  is  the  porosity  of  the 
sand,  C  [ML'  ]  is  the  concentration  of  dissolved  CAH  in  the  sand  aquifer  layer,  t  [T]  is 


time,  x  [L]  is  the  distance  along  the  respective  Cartesian  coordinate  axis,  Dfj  [L2!"1]  is  the 
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hydrodynamic  dispersion  coefficient  tensor  for  sand,  and  vf  [LT'1]  is  the  linear  pore 
water  velocity.  The  dispersion  coefficients  are  related  to  the  water  velocity  by  Equations 
3.2  and  3.3. 


Dl  =  aLvi 


(Eq  3.2) 


and 


Dx  =  axVi  (Eq3.3) 

where  Dl  [L2T_1]  is  the  longitudinal  dispersion  coefficient,  Dx  [L2T_1]  is  the  transverse 
dispersion  coefficient,  <Xl  [L]  is  the  longitudinal  dispersivity,  and  ax  is  the  transverse 
dispersivity.  Note  that  DL  is  Dfj  when  i  =  j  =  x,  and  x  is  the  direction  of  flow,  and  Dx  is 
Dfj  when  i  =  j  =  y  and  y  is  the  direction  transverse  to  the  flow  direction. 

In  Equation  3.1,  advection  in  the  sand  layer  is  represented  by  (0VjC),  and 
dispersion  is  represented  by  (  0Dfj  ^ ).  Equation  3.4  relates  the  sorbed  concentration 

0X[  y  •*  (7Xj  j 

to  the  concentration  of  dissolved  CAH  (C)  by  a  linear,  reversible,  equilibrium  isotherm: 

C  =  K^C  (Eq  3.4) 

where  C  [MM-1]  is  the  concentration  of  contaminant  sorbed  to  sand.  The  retardation 
factor  (R)  in  Equation  3. 1  is  defined  as: 

R=1  +  ff  =1+fKd  (Eq3.5> 

where,  p^fML"1]  is  the  bulk  density  of  sand,  and  [L’VT1]  is  the  sorption  distribution 
constant  for  sand. 

Darcy’s  Law,  represented  by  Equation  3.6  couples  the  hydraulic  conductivity  and 
gradient  to  the  advective  velocity  term. 
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qi  K{  Sh  .  , 

V‘=^  =  -J^  <Ec>3-6) 

where  q,  [L3L'2T‘1]  is  the  volumetric  flow  rate  per  unit  surface  area  of  aquifer  (also 
referred  to  as  the  Darcy  velocity),  K,  [LT1]  is  the  principal  component  of  the  sand 
hydraulic  conductivity  tensor,  and  h  [L]  is  the  hydraulic  head  in  the  sand  aquifer. 

Equation  3.7  is  Equation  3.1  with  the  addition  of  a  source/sink  term. 

R0^  =  A(eDs  _  A(ev  C) +  qsC  (Eq  3.7) 

at  dxi  V  1J  dxj  axj v  1  J  ~  v  M  ’ 

where  qs  [T1]  is  the  volumetric  flow  rate  per  unit  volume  of  aquifer  representing  fluid 

sources  (positive)  and  sinks  (negative). 

As  shown  in  Figure  4,  below  the  sand  aquifer  is  a  low  permeability  clay  layer  that 
contains  cracks  caused  by  pooled  DNAPL.  To  model  this  situation  a  dual  porosity,  or 
dual  domain,  model  is  used.  A  dual  domain  model  is  ideal  for  modeling  contaminant 
transport  in  fractured  media.  In  a  dual  domain  model,  water  is  assumed  to  be  either 
mobile  or  immobile.  Contaminant  can  be  transported  by  advection,  dispersion,  and 
sorption  in  the  mobile  water,  while  in  the  immobile  water,  there  is  no  contaminant 
transport,  though  equilibrium  sorption  may  be  modeled  using  an  immobile  domain 
retardation  factor.  Exchange  of  contaminant  between  the  mobile  and  the  immobile 
domains  is  assumed  to  be  governed  by  a  first-order  rate  process,  characterized  by  a  first- 
order  rate  constant  (Zheng  and  Wang,  1999). 

The  governing  equations  adapted  from  Zheng  and  Wang  (1999)  for  a  dual  domain 
system  with  sorption  are: 


0mRm 


dCn 

St 


+  SimRir 


SCi 


St 


Sxj  (0mDi:j  sx™)  Sxj  (0mViCCm)  (Eq  3.8) 


and 
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(Eq  3.9) 


0imRim^f  =  q(Cm-Cim) 

where  Cm  [ML'  ]  is  the  concentration  of  dissolved  CAH  in  the  mobile  domain, 

Cim  [ML  3]  is  the  concentration  of  dissolved  CAH  in  the  immobile  domain,  Rm  [-]  is  the 
retardation  factor  for  the  mobile  domain,  Rim  [-]  is  the  retardation  factor  for  the  immobile 
domain,  0m  [-]  is  the  porosity  of  the  mobile  domain,  defined  as  volume  of  mobile  water 
per  total  volume  of  the  medium,  0;m  [-]is  the  porosity  of  the  immobile  domain,  defined  as 
volume  of  immobile  water  per  total  volume  of  the  medium.  When  z=0  in  Figure  4  Cm 
and  Cim  equal  C  in  Equation  3.7.  Dfj  is  the  dispersion  coefficient  of  the  dissolved  CAH 
in  the  cracks,  vf  is  the  pore  water  velocity  in  the  cracks,  and  C  [T‘ 1  ]  is  the  first  order  mass 
transfer  rate  between  the  mobile  and  immobile  domains.  The  equation  for  Rm  is: 

Rm  =  l  +  hf^  (Eq  3.10) 

°m 

where  f  [-]  is  the  fraction  of  sorption  sites  that  are  in  contact  with  the  mobile  liquid  phase 
and  is  approximated  by  p£  is  the  bulk  density  of  the  cracked  clay,  and  [L  'M1]  is 
the  sorption  constant  for  clay.  The  equation  for  R;m  is: 

Rim  =  1  +  (1~eiP^  (Eq  3.11) 

uim 

The  source/sink  term,  which  is  the  last  tenn  on  the  right-hand  side  of  Equation  3.7,  is 
used  to  couple  transport  in  the  sand  with  Equations  3.8  and  3.9  describing  transport  in  the 
cracked  clay.  The  model  developed  for  this  research  assumes  that  no  biological  or 
chemical  degradation  of  the  contaminant  occurs. 
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3.4  Model  Implementation  in  GMS 

Using  DoD’s  Groundwater  Modeling  System  (GMS)  “enhanced  diffusion”  of 
dissolved  phase  CAH  into  low  penneability  clay  from  a  pool  of  DNAPL  and  back 
diffusion  of  dissolved  phase  CAH  out  of  the  low  penneability  clay  is  modeled.  GMS 
uses  MODFLOW  to  generate  the  model’s  groundwater  flow  contours.  MODFLOW 
determines  these  contours  by  first  generating  a  steady-state  potential  field  based  on 
boundary  conditions  and  the  main  equation  of  flow.  Then,  using  Equation  3.6,  with  the 
specified  hydraulic  conductivities,  MODFLOW  detennines  steady  groundwater  velocities 
throughout  the  model  domain.  The  velocity  field  is  then  used  in  Equation  3.7  or  3.8 
(depending  on  whether  a  single  or  dual  domain  transport  model  is  being  formulated). 
Equation  3.7  or  3.8  is  simulated  by  Modular  Three-Dimensional  Multispecies  Transport 
Model  Multi-Species  (MT3DMS)  to  simulate  advection/dispersion,  sorption,  diffusion 
and  biological  or  chemical  reactions  of  dissolved  contaminants  in  groundwater. 

Within  GMS,  a  100  m  long  by  70  m  wide  by  6  m  deep  sand  aquifer  was 
constructed  above  a  100  m  long  by  12  m  deep  by  30  m  wide  low  conductivity  clay  layer. 
In  MODFLOW  values  of  the  vertical  and  horizontal  hydraulic  conductivity,  longitudinal 
dispersivity  and  mobile  porosity  were  assigned  to  the  sand  and  clay.  Also  to  establish  a 
hydraulic  gradient  in  the  model,  the  heads  on  the  left  and  right  side  of  the  model  were 
20.5  m  and  19.5  m,  respectively.  Using  MT3DMS,  the  advection,  dispersion,  sorption, 
and  diffusion  transport  mechanisms  in  the  sand,  and  clay  are  simulated.  The  simulation 
consisted  of  two  stress  periods.  During  the  first  stress  period  of  10  years,  a  constant 
concentration  source  of  1 10  mg/L  of  TCE  was  introduced  just  above  a  192  m  area  of  the 
clay.  The  second  stress  period  was  20  years  long  during  which  no  additional  TCE  is 
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introduced.  MT3DMS  uses  a  third-order  TVD  (ULTIMATE)  scheme  to  solve  the 
advection  piece  of  the  model.  The  ratio  of  horizontal  transverse  dispersivity  to 
longitudinal  dispersivity  is  0.1  and  the  ratio  of  vertical  transverse  dispersivity  to 
longitudinal  dispersivity  is  0.01.  An  effective  molecular  diffusion  coefficient  of  8.64  x 
10‘6  L^T"1  is  assigned  to  the  sand  and  clay.  Dual-domain  mass  transfer  (with  sorption) 
was  assumed  for  every  model  run  except  Scenario  1  mentioned  below.  TCE 
concentration  versus  time  data  were  calculated  at  a  single  observation  point  56  m 
downgradient  from  the  center  of  the  TCE  source. 

Figure  5  is  a  conceptual  depiction  of  the  GMS  model.  The  non-aqueous  phase 
TCE  Source  (DNAPL  Pool)  is  transformed  into  aqueous  phase  TCE  and  is  transported 
into  the  mobile  domain  (crack)  and  downgradient  in  the  high  permeability  layer  (sand 
layer)  by  advection/dispersion.  The  TCE  source  is  represented  in  GMS  as  constant 
concentration  cells,  which  are  set  to  zero  when  the  source  is  removed.  Diffusion  between 
the  mobile  domain  and  the  immobile  domain  (low  permeability  clay)  is  modeled  as  a 
first-order  process  characterized  by  a  first  order  mass  transfer  rate  constant  (Q.  Within 
the  mobile  and  immobile  domains  the  aqueous  TCE  is  assumed  to  be  well  mixed.  Cm(z,t) 
and  Cim(z,t)  represent  the  aqueous  TCE  concentration  at  a  given  point,  z,  and  time,  t,  in 
the  mobile  and  immobile  domains,  respectively.  Both  the  mobile  and  immobile  domains 
have  retardation  factors  represented  by  Rm  for  the  mobile  domain  and  R;m  for  the 
immobile  domain.  Rm  and  R,m  are  calculated  using  Equations  3.10  and  3.11  respectively. 
Ed  represents  the  ratio  of  sorbed  to  dissolved  TCE  concentrations  in  the  clay. 
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Figure  5:  Conceptual  Depiction  of  Dual-Domain  GMS  Model. 


To  simulate  cracks  in  the  clay,  the  values  of  the  clay’s  mobile  and  immobile 
porosity  were  modified  in  MODFLOW  and  MT3DMS.  To  simulate  the  formation  of 
additional  cracks  the  clay’s  mobile  porosity  (0m)  was  increased  as  the  immobile  porosity 
(0jm)  decreased.  In  the  absence  of  cracks,  the  clay  was  assumed  to  have  a  mobile  porosity 
of  zero  and  an  immobile  porosity  (0^)  of  0.43.  To  calculate  the  immobile  porosity  (0im) 
of  the  cracked  clay,  Equation  3.12  is  used. 

Qim  =  e?m(l  -  9m)  (Eq3.12) 

where  0m  [-]  is  the  mobile  porosity  assigned  in  MODFLOW  and  0;m  [-]  is  assigned  in 
MT3DMS.  To  determine  the  total  porosity  (0)  of  cracked  clay  Equation  3.13  was  used. 

0  =  0m  +  OfmCl  -  0m)  =  em  +  eim  (Eq  3.13) 

In  MT3DMS  values  of  bulk  density  values  were  assigned  to  the  sand  and  clay. 

The  bulk  density  (pb)  was  calculated  using  Equation  3.14. 
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Pb  =  Ps(l  -  6) 


(Eq  3.14) 


3  •  •  •  •  3 

where  ps  [kg-trf  ]  is  the  solid  density.  The  value  of  ps  for  sand  is  2649  kg-nT  and  for 
clay  is  2630  kg-nT  (Chapman  and  Parker,  2005).  The  0  of  the  sand,  which  is  set  equal  to 
0m  of  the  sand  since  sand  is  assumed  to  have  no  immobile  porosity,  is  0.35  (Chapman  and 
Parker,  2005).  Thus,  the  sand  pb  value  used  in  every  model  run  is  1722  kg-nC .  The 
value  of  pb  assigned  to  the  clay  varied  from  model  run  to  model  run  because  the  value  of 
Pb  depends  on  the  value  of  0,  per  Equation  3.14,  and  the  value  of  0  varies  with  the  value 
of  0m  assigned  to  the  clay,  per  Equation  3.13. 

To  solve  for  B  in  Figure  4  Equation  3. 15  is  used  (Goltz  and  Roberts,  1988). 


?  _  3Deff 


0imRim  B2Rim 


(Eq  3.15) 


where  B  [L]  is  the  half-distance  between  cracks  as  depicted  in  Figure  4,  ^  [T1]  is  the  rate 
constant  for  mass  transfer  between  the  low  permeability  clay  and  the  crack,  0im  [-]  is  the 
porosity  of  the  low  permeability  clay,  Rim  [-]  is  the  retardation  factor  of  the  low 
permeability  clay,  and  Deff  [L2T'1]  is  the  effective  diffusion  coefficient  within  the  clay 
calculated  using  Equation  2.2. 

To  solve  for  the  radius,  b,  of  a  single  crack  the  total  volume  of  cracks  must  first 
be  detennined.  Equation  3.16  is  used  to  detennine  the  total  volume  of  cracks  which  is 
represented  by  the  volume  of  mobile  void  space. 


0 


m 


volume  mobile  void  space 
Total  Volume 


(Eq  3.16) 


where  0ra  [-]  is  the  porosity  of  the  mobile  void  space  in  the  low  penneability  layer.  The 
cracks  in  the  low  penneability  clay  represent  the  mobile  void  space.  With  the  total 
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volume  of  cracks  determined,  Equation  3.17  can  be  used  to  determine  the  radius,  b  of  a 
single  crack. 


Volume  Mobile  Void  Space  =  CTib2h  (Eq  3.17) 

where  c  [-]  is  the  number  of  cracks  in  the  clay,  b  [L]  is  the  radius  of  a  single  crack,  and  h 
[L]  is  the  height  of  the  crack  (see  Figure  4).  The  number  of  cracks  in  the  clay,  needed  to 
apply  Equation  3.17,  is  determined  using  Equation  3.18. 

c=(jj+1)(j5+!)  (Eq  3.18) 

where  L  [L]  is  the  length  and  W  [L]  is  the  width  of  the  cracked  clay  region. 

Table  1  lists  the  crack  radius  (b)  as  depicted  in  Figure  4,  and  the  half  distance  between 
cracks  (B)  calculated  using  Equations  3.15  through  3.18  based  on  the  value  of  0m. 


Table  1:  Crack  Radius,  Half  Distance  Between  Cracks  and  Total  Number  Cracks  Based  On 


6m  (") 

c(-) 

b(m) 

B  (m) 

0.001 

0 

0 

0 

0.01 

341,239 

0.005 

0.047 

0.1 

375,292 

0.016 

0.045 

0.5 

674,859 

0.027 

0.033 

3.5  Model  Scenarios 

Two  initial  model  scenarios  were  considered.  Scenario  1  was  a  high  permeability 
sand  that  contained  a  TCE  source.  Table  2  lists  the  MODFLOW  parameters  used  and 
Table  3  lists  the  MT3DMS  parameters  used.  Due  to  the  fact  that  there  was  no  clay  in  this 
scenario  the  model  was  run  without  immobile  porosity  or  first  order  mass  transfer  rate 
constant. 
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Table  2:  Scenario  1  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

N/A 

N/A 

N/A 

N/A 

Ul  (Chapman  and 

.  Parker,  2005) 

Table  3:  Scenario  1  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(pb)  (kg-nT3) 

Immobile  Porosity 

(0im) 

Sorption  Constant 
(Kd)  (m3-kg_1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

N/A 

0.00017(1) 

N/A 

Clay 

N/A 

N/A 

N/A 

N/A 

(1)(Sale  et  al,  2008) 


Scenario  2  was  high  permeability  sand  located  above  a  low  permeability  clay 
layer  that  was  free  of  cracks.  On  the  surface  of  the  low  permeability  clay  layer  was  a 
pool  of  TCE.  Table  4  lists  the  MODFLOW  parameters  used  and  Table  5  lists  the 
MT3DMS  parameters  used.  As  mentioned  earlier  this  research  assumed  that  uncracked 
clay  has  no  mobile  porosity.  To  simulate  this,  a  value  of  0.001  was  assigned  for  the 
mobile  porosity  of  the  clay  (see  Table  4).  Recall  the  clay  total  porosity  is  0.43;  thus,  an 
immobile  porosity  of  0.429  was  assumed  for  the  uncracked  clay  (see  Table  5). 


Table  4:  Scenario  2  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.001 

(1)  (Chapman  and 

.  Parker,  2005) 
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Table  5:  Scenario  2  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(pb)  (kg-nT3) 

Immobile  Porosity 

(Oim) 

Sorption  Constant 
(Kd)  (m3-kg_1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1499 

0.429 

0.051 

0.005 

“’(Sale  et  al,  2008) 


The  concentration  versus  time  breakthrough  curves  generated  for  Scenarios  1  and 
2  are  used  as  “benchmarks”  for  the  parameter  sensitivity  analysis  discussed  later. 

Scenario  3  was  the  same  as  Scenario  2  except  Scenario  3  made  the  assumption 
that  vertical  cracks  existed  in  the  clay.  Values  of  IQ,  C,  and  0m  used  in  Scenario  3  are 
baseline  values  for  the  sensitivity  analysis  to  be  discussed  later.  Table  6  lists  the 
MODFLOW  parameter  values  used  and  Table  7  lists  the  MT3DMS  parameter  values 
used. 


Table  6:  Scenario  3  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.1 

(1)  (Chapman  and 

.  Parker,  2005) 

Table  7:  Scenario  3  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(Pb)  (kg-nT3) 

Immobile  Porosity 

(0im) 

Sorption  Constant 
(K„)  (m3-kg_1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1349 

0.387 

0.051 

0.005 

“’(Sale  et  al.,  2008) 
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3.6  Sensitivity  Analysis 


A  sensitivity  analysis  was  conducted  to  investigate  how  varying  the  sorption 
distribution  coefficient  (IQ),  mobile  porosity  (0m),  and  the  first  order  mass  transfer  rate 
constant  (Q  affected  the  downgradient  concentration  versus  time  breakthrough  curves. 
The  sensitivity  analysis  conducted  for  IQ  was  performed  done  using  values  of  0.005 1  m3- 
kg'1,  and  0.51  nr-kg'1.  The  breakthrough  curves  generated  are  compared  to  the 
breakthrough  curves  for  Scenarios  1-3.  Table  8  lists  the  MODFLOW  parameters  used 
and  Table  9  lists  the  MT3DMS  parameters  used. 


Table  8:  Kd  Sensitivity  Analysis  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.1 

Ul  (Chapman  and 

.  Parker,  2005) 

Table  9:  Kd  Sensitivity  Analysis  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(pb)  (kg-nT3) 

Immobile  Porosity 

(0im) 

Sorption  Constant 

(Kd)  (m3-kg1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1349 

0.387 

0.0051 

0.005 

0.51 

“’(Sale  et  al„  2008) 


The  sensitivity  analysis  for  0m  was  done  using  0m  values  of  0.01,  and  0.5.  The 
breakthrough  curves  generated  were  compared  to  the  breakthrough  curves  for  Scenarios 
1-3.  Table  10  lists  the  MODFLOW  parameter  values  used  and  Table  1 1  lists  the 
MT3DMS  values  used.  For  each  0m  value  in  Table  9  Equation  13.11  was  used  to 


42 


calculate  0im  and  Equation  13.13  was  used  to  calculate  the  corresponding  pb  value  in  Table 

11. 


Table  10:  0m  Sensitivity  Analysis  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

0.01 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.1 

0.5 

Ul  (Chapman  and 

.  Parker,  2005) 

Table  11:  8m  Sensitivity  Analysis  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(pb)  (kg-m3) 

Immobile  Porosity 
(Bin,) 

Sorption  Constant 
(Kd)  (nf-kg1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1484 

0.4257 

0.051 

0.005 

750 

0.215 

“’(Sale  et  al„  2008) 


The  final  sensitivity  analysis  conducted  was  for  C,  values  of  0.01  d'1,  and  0.5  d"1. 
The  breakthrough  curves  generated  are  compared  to  the  breakthrough  curve  for  Scenarios 
1-3.  Table  12  lists  the  MODFLOW  parameters  used  and  Table  13  lists  the  MT3DMS 
parameters  used. 


Table  12:  C,  Sensitivity  Analysis  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.1 

Ul  (Chapman  and 

.  Parker,  2005) 
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Table  13:  C,  Sensitivity  Analysis  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(pb)  (kg-m'3) 

Immobile  Porosity 

(0im) 

Sorption  Constant 
(Kd)  (m3-kg_1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1349 

0.387 

0.051 

0.01 

0.5 

(1)(Sale  et  al„  2008) 


3.7  First  Moment  Analysis 

Since  back  diffusion  causes  breakthrough  curve  tailing,  and  an  understanding  of 
the  extent  of  tailing  is  crucial  in  managing  a  contaminant  plume,  an  analysis  of  the 
breakthrough  curve’s  first  moment  was  conducted  for  this  study.  The  first  moment 
measures  the  center  of  mass  of  the  breakthrough  curve.  A  high  value  for  the  first 
moment  indicates  a  high  degree  of  tailing.  In  this  study,  the  first  moment  of 
breakthrough  curves  was  determined  for  IQ  values  of  0.051  m  -kg'  and  0.102  m  -kg'  . 
Table  14  lists  the  MODFLOW  parameters  used  and  Table  15  lists  the  MT3DMS 
parameters  used.  The  total  model  run  time  was  695  years  to  accurately  capture  the  first 
moment  of  the  tailing  breakthrough  curves.  The  constant  1 10  mg-L'1  TCE  source  was 
removed  after  10  years. 


Table  14:  Kd  Moment  Analysis  MODFLOW  Material  Property  Input  Parameters 


Material 

Horizontal  Kh 
(m/d) 

Vertical  Kv 
(m/d) 

Longitudinal 
Dispersivity  (m) 

Mobile 
Porosity  (0m) 

Sand(1) 

17.28 

1.728 

1.0 

0.35 

Clay 

4.32E-5(1) 

4.32E-6(1) 

1.0E-4 

0.1 

(1)  (Chapman  and 

.  Parker,  2005) 
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Table  15:  Kd  Moment  Analysis  MT3DMS  Input  Parameters 


Material 

Bulk  Density 
(Pb)  (kg-nT3) 

Immobile  Porosity 
(Oim) 

Sorption  Constant 
(Kd)  (m3-kg_1) 

First  Order 
Mass 

Transfer  Rate 
Constant  (Q 
(d1) 

Sand 

1722 

0.00001 

0.00017(1) 

0.0 

Clay 

1349 

0.387 

0.051 

0.102 

0.005 

'(Sale  et  al„  2008) 


To  calculate  the  first  moment,  the  area  under  the  concentration  versus  time  curve 
was  calculated  using  Equation  3.18. 

area  =  [c(ti+1  J+C(ti)]  (tj+1  _  tj)  (Eq  3.18) 

where  c(t)  [mg-L'1]  is  the  concentration  output  from  GMS  and  tj  [T]  is  the  time. 

Then  the  discrete  residence  time  density  function  (f(h))  is  determined  for  each  discrete 
concentration  value  using  Equation  3.19. 

f(ti)=^  (Eq  3.19) 

The  first  moment  is  then  calculated  using  Equation  3.20. 

T  _  V*max— 1  [(h+h+i)l  [hfO+hfi+i)!  r*-  +-  \ 

tRTD  —  Li=o  2  2  ]  V-i+l  m)  (Eqi./U) 

where  tRTD  [T]  is  the  mean  residence  time  or  first  moment  (Clark,  2009). 
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4.0  Results  and  Discussion 


4.1  Overview 

This  chapter  presents  the  breakthrough  curves  for  Scenarios  1-3  and  the  IQ,  0m, 
and  £  sensitivity  analysis  discussed  in  Chapter  3.  The  concentration  versus  time  data  was 
calculated  at  a  single  observation  point  56  m  downgradient  from  the  TCE  source.  Also 
presented  are  the  results  of  the  first  moment  analysis  conducted  for  IQ  values  of  0.05 1 
nr -kg'1  and  0.102  m3-kg"' . 

4.2  Baseline  Simulations 

Figure  6  is  the  breakthrough  curves  for  baseline  Scenarios  l-3along  with  the  TCE 
maximum  contaminant  level  (MCL)  of  0.005  mg-L'1  for  comparison. 
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Figure  6:  Breakthrough  Curves  for  Scenarios  1-3 
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These  three  breakthrough  curves  will  be  used  to  compare  to  the  breakthrough  curves  for 
the  Kd,  0m,  and  C  sensitivity  analysis. 

Scenario  1  represents  high  penneability  sand  with  a  constant  110  mg-L'1  TCE 
source  that  disappears  after  10  years.  For  the  remaining  20  years,  the  sand  is  flushed  with 
uncontaminated  water.  The  breakthrough  curve  in  Figure  6  indicates  that  after  the  TCE 
source  is  removed  the  aqueous  phase  TCE  concentration  in  the  sand  rapidly  drops  and 
tails  off  towards  zero.  The  tailing  effect  is  not  represented  in  Figure  6. 

Scenario  2  consists  of  high  permeability  sand  atop  an  uncracked  low  permeability 
clay  layer.  A  1 10  mg-L'1  TCE  source  is  used  to  represent  a  DNAPL  pool  atop  the  clay 
layer.  The  source  remains  for  10  years.  This  scenario  assumes  that  the  dominant 
transport  mechanism  of  the  aqueous  phase  TCE  into  the  clay  from  the  DNAPL  pool  is 
Fickian  diffusion.  The  tailing  of  the  breakthrough  curve  beginning  at  approximately  1 1 
years  for  Scenario  2  reflects  the  “back  diffusion”  of  the  aqueous  phase  TCE  from  the 
uncracked  clay  into  the  sand  aquifer  where  advection  then  becomes  the  dominant 
transport  mechanism.  Figure  6  also  indicates  that  the  downgradient  aqueous  phase  TCE 
concentration  is  still  above  the  U.S.  EPA  maximum  contaminant  level  (MCL)  of  0.005 
mg-L'1  20  years  after  the  TCE  source  is  removed. 

Scenario  3  is  the  same  as  Scenario  2  except  Scenario  3  contains  cracked  clay 
instead  of  just  clay.  The  tailing  in  Scenario  3  is  greater  than  the  tailing  in  Scenario  2,  due 
to  enhanced  “back  diffusion”  of  aqueous  phase  TCE  from  the  cracked  clay  into  the  sand 
aquifer.  Also,  the  downgradient  TCE  concentration  is  still  above  the  U.S.  EPA 
maximum  contaminant  level  (MCL)  of  0.005  mg-L'1  20  years  after  the  TCE  source  is 
removed. 
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Tailing,  which  is  the  significant  difference  between  the  three  breakthrough  curves 
in  Scenarios  1  through  3,  is  due  to  the  low  permeability  clay  in  Scenarios  2  and  3.  The 
aqueous  phase  TCE  is  “stored”  in  the  clay  until  the  DNAPL  source  is  depleted  in  the  sand 
aquifer.  Upon  removal  of  the  DNAPL  source,  the  aqueous  phase  TCE  “stored”  in  the 
clay  “back  diffuses”  into  the  sand,  where  it  is  then  transported  downgradient  by 
advection.  This  “back  diffusion”  process  results  in  the  breakthrough  curve  tail  for 
Scenarios  2  and  3.  The  tailing  that  occurs  in  Scenarios  2  and  3  is  similar  to  the  tailing 
observed  by  Chapman  and  Parker  (2005)  in  their  investigation  of  back  diffusion  of  TCE 
from  a  clayey  silt  aquitard. 

4.3  Sensitivity  Analysis 

A  sensitivity  analysis  is  conducted  for  Kd,  0m,  and  C.  For  each  sensitivity  analysis, 
the  breakthrough  curve  is  plotted  and  compared  to  the  breakthrough  curves  generated  for 
the  baseline  scenarios.  The  first  sensitivity  analysis  conducted  considers  how  the  results 
of  the  model  would  change  if  Kd  of  the  clay  is  changed.  The  values  of  Kd  considered  are 
0.005 1  m  -kg'  ,  and  0.5 1  m  -kg'  .  Baseline  parameter  values  for  9m  and  C  oft).  1  and 
0.005  d'1  respectively  were  used.  Figure  7  depicts  the  breakthrough  curves  for  the  two 
different  Kd  values  along  with  the  breakthrough  curves  the  baseline  scenarios. 
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Figure  7:  Breakthrough  Curves  for  Kd  Sensitivity  Analysis  along  with  Baseline  Scenarios 


The  distribution  coefficient,  IQ,  represents  the  partitioning  of  the  contaminant 
between  the  solid  and  the  aqueous  phase.  As  the  value  of  IQ  increases  the  concentration 
of  the  contaminant  in  the  sorbed  phase  increases  and  the  concentration  of  the  contaminant 
in  the  aqueous  phase  decreases.  In  the  model,  as  the  value  of  the  clay  Kd  is  increased  the 
concentration  of  the  dissolved  contaminant  measured  at  the  downgradient  observation 
point  decreases  because  more  of  the  contaminant  is  in  the  sorbed  phase  in  the  low 
permeability  clay  layer.  This  increase  in  sorbed  phase  concentration  leads  to  a  decrease 
in  the  “back  diffusion”  of  the  aqueous  phase  out  of  the  low  permeability  clay  after  the 
original  TCE  source  is  depleted.  This  can  be  seen  in  the  breakthrough  curves  plotted  in 
Figure  7  for  the  three  values  of  Kd  chosen;  as  Kd  increases,  the  concentration  at  which 
“tailing”  occurs  decreases.  In  fact,  for  the  highest  value  of  IQ  used  in  the  analysis,  the 
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concentration  of  the  tail  is  below  the  MCL  for  TCE.  Note,  however,  that  since  more 
contaminant  is  stored  in  the  low  conductivity  layer  when  Kd  is  increased,  ultimately  more 
will  be  released.  This  effect  will  be  quantified  when  we  calculate  the  first  moment  of  the 
breakthrough  curves  for  different  values  of  Kd  in  Section  4.4. 

The  next  sensitivity  analysis  examines  the  first  order  mass  transfer  rate  constant 
(Q.  Values  of  0.01  d'1,  and  0.5  d'1  are  chosen  and  the  breakthrough  curves  generated  are 
compared  to  the  breakthrough  curves  for  the  baseline  scenarios.  Baseline  parameter 
values  for  Kd  and  0m  of  0.05 1  nr  -kg'  and  0. 1  respectively  were  used.  Figure  8  depicts 
the  breakthrough  curves  for  the  C,  sensitivity  analysis  along  with  the  breakthrough  curves 
for  the  baseline  scenarios. 
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Figure  8:  Breakthrough  Curves  for  C,  Sensitivity  Analysis  along  with  Baseline  Scenarios 

(MT  =  -Mass  Transfer) 
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C,  represents  the  mass  transfer  rate  between  the  mobile  and  immobile  domains.  As 
the  value  of  C,  is  increased,  contaminant  mass  is  transferred  more  quickly  from  the  cracks 
in  the  clay  (the  mobile  domain)  to  the  clay  (the  immobile  domain).  This  increase  in  mass 
transfer  leads  to  an  increase  in  contaminant  mass  stored  in  the  clay  which  results  in  a 
lower  aqueous  phase  TCE  concentration  downgradient  after  the  TCE  source  is  depleted. 
As  C,  decreases  and  mass  transfer  slows,  more  mass  remains  in  the  mobile  domain,  so 
higher  concentrations  of  dissolved  TCE  back  diffuse  initially,  when  the  DNAPL  source  is 
removed.  Ultimately,  as  may  be  seen  qualitatively  by  looking  at  Figure  8,  since  the  value 
of  Ed  is  the  same  for  the  three  values  of  the  first  moments  of  the  breakthrough  curves 
are  the  same.  This  phenomenon  has  been  recognized  by  a  number  of  investigators  (e.g., 
Valocchi,  1985;  Goltz  and  Roberts,  1986). 

The  last  sensitivity  analysis  varies  the  mobile  porosity,  0m,  of  the  clay.  Mobile 
porosity  values  of  0.01,  and  0.5  are  chosen  and  the  breakthrough  curves  generated  are 
compared  to  the  breakthrough  curves  for  the  baseline  scenarios.  Baseline  parameter 
values  for  IQ  and  and  C,  of  0.05 1  m  -kg'  and  0.005  d'  respectively  were  used.  Figure  9 
depicts  the  breakthrough  curves  for  the  0m  sensitivity  analysis. 


51 


0) 

o 

s 

s 

U 

U 


1.0E-01 


1.0E-02 


1.0E-03 


1.0E-04 


1.0E-05 


- Tketa_m=0.01 

- Theta  m=0.5 

- Scenario  1-  Sand  Only 

. Scenario  2-  Sand  and  Clay 

. Scenario  3-  Sand  and  Cracked  Clay 

-  TCE  MCL  (0.005  mg/L) 

\ 

i  V. 

- , — 

0  5  10  15  20  25  30  35 


Time(y) 


Figure  9:  Breakthrough  Curves  for  Clay  0m  Sensitivity  Analysis  along  with  Baseline 

Scenarios 


For  this  research  cracks  in  the  clay  are  represented  by  0m.  As  the  value  of  the  clay 
0m  increases  the  number  of  cracks  in  the  clay  increases.  The  cracks  in  the  clay  lead  to 
“enhanced  diffusion”  of  the  aqueous  phase  TCE  into  the  clay  from  the  DNAPL  source  on 
the  surface  of  the  clay.  This  “enhanced  diffusion”  process  ultimately  leads  to  an  increase 
in  the  contaminant  mass  stored  in  the  clay  as  compared  to  the  contaminant  mass  stored  in 
uncracked  clay.  Upon  depletion  of  the  DNAPL  source,  the  contaminant  mass  is 
transported  from  the  clay  into  the  sand  aquifer  by  “back  diffusion”  and  then 
downgradient  by  advection.  This  is  demonstrated  by  comparing  the  breakthrough  curve 
for  Scenario  2  which  is  the  breakthrough  curve  for  uncracked  clay  (0m  =  0.001)  and  the 
breakthrough  curve  for  clay  with  substantial  cracking  (0m  =  0.5).  The  tail  for  the 
breakthrough  curve  when  0m  equals  0.5  is  at  a  greater  concentration  than  the  tail  of  the 
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breakthrough  curve  for  the  uncracked  clay.  An  increase  in  9m  simulates  an  increasing 
number  of  cracks,  and  an  increasing  number  of  cracks  results  in  an  increase  in 
contaminant  mass  transported  into  and  stored  in  the  clay  (i.e.,  enhanced  diffusion).  In 
turn,  the  enhanced  diffusion  leads  to  increased  tailing  as  contaminant  that  back  diffuses 
from  the  clay  into  the  sand  aquifer  upon  depletion  of  the  DNAPL  source  is  transported 
downgradient  by  advection. 

Table  1  in  Chapter  3  provides  supporting  evidence  for  the  idea  that  as  0ra 
increased  the  number  of  cracks  in  the  clay  increased.  As  Table  1  shows  as  0m  increased 
from  0.01  to  0.5  the  number  of  cracks  in  the  clay  increased  from  341,239  to  375,292. 

Also  shown  in  Table  1  is  the  fact  that  as  0m  increased  from  0.01  to  0.5  the  radius  of  the 
cracks  (b)  increased  from  0.005  m  to  0.027  m  and  the  half  distance  between  the  cracks 
(B)  decreased  from  0.047  m  to  0.033  m.  These  results  are  consistent  with  the  fact  that  as 
9m  of  the  clay  increased  more  cracks  formed  and  the  half  distance  between  the  cracks 
must  decrease  to  accommodate  the  increasing  number  of  cracks  formed. 

4.4  First  Moment  Analysis 

A  moment  analysis  was  conducted  and  the  results  compared  for  Kd  values  of 
0.051  m  -kg'  and  0.102  m  -kg'  .  The  analysis  used  the  baseline  values  for  0ra  and  C,  of 
0.1  and  0.005  d'1  respectively.  This  analysis  was  conducted  to  verify  that,  the  appearance 
of  Figure  7  notwithstanding,  ultimately  the  extent  of  tailing  of  a  breakthrough  curve  with 
a  high  Kd  is  greater  than  for  a  lower  Kd.  Mathematically  this  may  be  accomplished  by 
calculating  and  comparing  the  values  of  the  first  moment  for  each  Kd  value. 
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The  first  moment  represents  the  mean  residence  time  of  the  solute  molecules  in 
the  system.  As  the  value  of  Kd  increases  more  solute  is  in  the  sorbed  phase  and  less 
solute  is  in  the  dissolved  aqueous  phase.  This  leads  to  an  increased  solute  residence  time 
in  the  system  and  a  higher  first  moment  value.  As  the  value  of  Kd  decreases  more  solute 
is  in  the  dissolved  aqueous  phase  which  leads  to  a  decrease  of  the  residence  time  in  the 
system  and  a  lower  first  moment  value. 

As  expected,  and  belying  the  appearance  of  Figure  7,  for  a  Kd  of  0.05 1  m  -kg'  the 
first  moment  is  approximately  1,113  years  compared  to  the  first  moment  of  1,259  years 
for  a  Kd  of  0.102  rrr’-kg"1. 
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5.0  Conclusion  and  Recommendations 


5.1  Conclusion 

The  primary  objective  of  the  research  was  to  model  the  effects  that  cracks  have  on 
the  storage  and  transport  of  dissolved  aqueous  phase  CAHs  into  and  out  of  the  low 
permeability  layers,  and  simulate  how  this  might  affect  the  evolution  and  longevity  of 
dissolved  phase  plumes.  The  dual-domain  model  developed  in  this  work  was  used  to 
simulate  cracks  that  might  form  in  the  low  permeability  clay  due  to  the  presence  of 
DNAPL.  These  cracks  allow  “enhanced  diffusion”  of  aqueous  phase  CAHs  into  the  clay 
which  leads  to  increased  contaminant  mass  stored  in  the  clay.  This  mass  stored  in  the 
clay  then  “back  diffuses”  out  of  the  clay  after  the  DNAPL  source  is  removed  creating  a 
dissolved  phase  contaminant  plume  that  exists  long  after  the  DNAPL  source  is  removed. 
Based  on  the  model  results,  we  found  that: 

1)  as  the  number  of  cracks  in  the  clay  increases,  “back  diffusion”  out  of  the  clay  after 
the  DNAPL  source  is  removed  increases,  and  a  higher  downgradient  concentration 
tail  persists  long  after  the  DNAPL  source  is  removed. 

2)  as  the  aqueous  phase  concentration  in  the  clay  increases  and  the  sorbed  phase 
concentration  decreases  a  higher  downgradient  concentration  tail  persists  long  after 
the  DNAPL  source  is  removed  due  to  increased  “back  diffusion”  out  of  the  clay. 

3)  when  the  sorbed  phase  concentration  in  the  clay  is  significantly  higher  than  the 
aqueous  phase  concentration,  “back  diffusion”  out  of  the  clay  is  not  significant  which 
could  lead  to  measured  downgradient  concentrations  that  are  below  regulatory  limits. 
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4)  upon  removal  of  a  DNAPL  source  pooled  atop  a  low  penneability  clay  layer  the 
“back  diffusion”  of  the  aqueous  phase  CAH  out  of  the  clay  can  cause  downgradient 
plume  persistence  long  after  the  DNAPL  source  is  removed. 

5)  remediation  strategies  that  focus  only  on  the  removal  of  the  DNAPL  source  and 
do  not  take  into  account  “back  diffusion”  out  of  the  low  permeability  clay  may  fail  to 
meet  site  cleanup  targets  and  goals. 

6)  the  model  has  the  potential  to  be  used  by  site  remediation  managers  to  guide  them 
in  their  decision  making  process  on  ways  to  best  manage  sites  requiring  remediation. 

5.2  Recommendations  for  Future  Research 

The  current  research  assumes  that  aqueous  phase  CAH  is  transported  from  the 
DNAPL  source  into  the  cracks  in  the  clay  by  advection.  Future  work  should  consider  the 
possible  impact  of  DNAPLtransport  into  the  cracks.  Clearly,  this  will  result  in  more 
contaminant  mass  in  the  low  penneability  layers.  The  plausibility  of  this  scenario 
depends  upon  such  factors  as  the  entry  pressure  required  for  DNAPL  to  enter  a  crack, 
which  in  turn  would  depend  upon  the  crack  diameter. 

Another  recommendation  is  to  use  mobile  porosity  and  hydraulic  conductivity 
values  for  clay  from  an  actual  site  that  has  been  in  contact  with  “real”  DNAPL 
(containing  impurities).  These  experimental  values  could  be  used  as  input  parameters  in 
the  model  developed  for  this  study. 

Other  recommendations  include  incorporation  of  contaminant  degradation  in  the 
clay  which  could  lead  to  a  lower  contaminant  concentration  that  “back  diffuses”  out  of 
the  clay,  two-dimensional  laboratory  experiments  to  validate  the  model,  investigation  of 
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other  models  that  could  be  used  to  simulate  cracking  in  the  clay,  development  of  a  model 
that  would  allow  for  cracks  to  fonn  over  time,  rather  than  assume  that  the  cracks  develop 


instantaneously  upon  arrival  of  the  contaminant  source,  and  finally  modeling  a  “real 
world”  remediation  scenario  would  help  with  model  validation. 
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